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Abstract. The Spitzer-GLIMPSE point source catalog and images have been 
used to study a sample of 381 massive protostellar candidates. IRAC-Point 
source photometry was used to analyse colours, magnitudes and spectral indicies 
of the infrared counterparts (IRCs) of massive protostellar objects and a bonafide 
sample of 50 point sources was obtained. Spectral energy distributions (SEDs) 
of these 50 sources was extended to the near-infrared and millimeter range by 
using 2MASS and millimeter data from the literature. An online SED fitter 
tool based on Monte-Carlo radiative transfer of an accretion model involving 
star, disk and envelope was used to fit the SEDs of the 50 sources. The IRCs 
to massive protostellar objects are found to successfully imitate the SEDs of 
evolutionary phases similar to low mass star formation. Envelope accretion, 
rather than disk accretion is found to be dominant in building the most massive 
stars. Unresolved centimeter continuum emission is associated with 27 IRCs 
classified as massive protostars suggesting that ionised accretion flows may play 
an important role along with the molecular component. The morphology of 
the infrared nebulae surrounding the IRCs have an unusual resemblance to the 
morphologies of ultra-compact HII regions suggesting that these infrared nebulae 
are possible precursors to the UCHII regions. 



1. Introduction 

In recent years, extensive observational data in the infrared and millimeter bands 
have led to the belief that the formation of massive stars mostl y involves accre - 
tion phenomenon similar to the formation of low mass stars (e.g. lWhitnev 2005). 



Whether the massive protostars are scaled up versions of low mass protostars 
or if they go through continuuing accretion on t o zero-age main sequence stars 
is an issue under debate (e.g: see the review bv lZinnecker & Yorkd 120071 ') . Ad- 



dressing these issues through direct observations requires angular resolutions of 
an order of magnitude better than what is currently possible in the infrared and 
millimeter range. The recently made available GLIMPSE survey already pro- 
vides a breakthrough in the spatial scales and sensitivity which is much higher 
than the upto now popular MSX and IRAS surveys. The Spitzer-IRAC point 
spread function of ~2" in the 3.6-8 fim bands are comparable to that of 2MASS 
in the J,H,K bands which allows us to probe scales of a few thousand AU at 
mean distances of 3-5 kpc where high mass protostellar objects (HMPOs) are 
commonly found. Similarly 1" resolution observations in the millimeter bands 
made possible by the Plateau de Bure interferometer or the Sub Millimeter Ar- 



ray has unveiled disks an d toroids around several HMPOs (e.g. iBeltran et al 
20061 : iBeuther et all 120071 ) . 
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A decade long effort, that used systematic selection criteria involving con- 
straints on the IRAS colours and other signposts of massive star formation such 
as masers and outflows, has resulted in four HMPO can didate samples, cov- 
ering the northern and southern hemispheres of the sky dMolinari et al" 19961 : 
Sridharan et~aDl2002l : iFontani et aUl2002l : iFaundez et alJl2004l ) (hereafter Mol96. 



Sri02, Fon02, Faun04). These surveys have provided a sample of ~ 500 objects 
and a major fraction of the sample has been investigated extensively to evaluate 
the dense gas, dust continuum, ionis ed emission and ma ser associations. Us- 
ing the 2MASS point source catalog, iKumar et al.l (2006) looked for clustering 
around the HMPOs in the northern hemisphere and found 54 embedded clusters 
and several near-infrared counterparts to these objects. In this work, we extend 
such as study to a more complete sample of HMPOs using the GLIMPSE survey 
data. 



2. Infrared counterparts to HMPOs 

Of the 500 targets from the Mol96, Sri02, Fon02 and Faun04, 381 targets were 
covered by the GLIMPSE fields for which point source photometry and images 
were extracted. Point sources within a 100" radius centered on each target 
were used for colour, magnitudes and spectral index analysis whereas images 
with a size of 300" were extracted for examining nebulosities and clustering. 
The distance to each target was taken as listed in the original references and 
a mean distance of 5kpc was assumed for sources when an estimate was not 
available. The point sources were placed on a [3.6-4.5] vs [5.6-8.0] colour-colour 
diagram and separat ed into the evolu t ionar y stages of Class I, II and III based 
on the definitions of iRobitaille et al.l (j2006j) (hereafter RWIWD06). The fluxes 



in the four IRAC bands namely 3.6, 4.5, 5.8 and 8.0 /mi bands were fitted with a 
simple least squares linear fit to obtain the spectral index a in these bands. Since 
the interest is not only to finding redenned objects but also luminous objects, we 
defined a parameter called alpha- magnitude product AM = -Mg Mm x(a-|-6)/10 
where Mg^ m is the 8.0/xm absolute magnitude of the source and a the observed 
IRAC spectral index. The constants 6 and 10 in the above equation were chosen 
arbitrarily to effectively separate the high a sources from the rest. For more 
details see IKumar fc Gravel (|2007n . 



In Fig. la we show a [3.6-4.5] vs [5.6-8.0] colour-colour diagram for all the 
sources detected from all target lists. The main concentration of points at (0,0) 
represents photospheres and higher values on both the x and y axes represent 
young stellar objects (YSOs) considered to be at different evolutionary stages. 
The regions representative of the evolutionary stages such as I, and II based on 
the 2D radiative transfer model data of RWIWD06 are marked. The diamonds 
represent candidate massive protostars which are those points lying to the right 
side of the vertical dashed line in Figlb. Fig. lb shows a AM product vs Mg^ m 
(absoulte magnitude) for point sources from all the targets. The solid curves 
represent 20M Q and 8M class I object models for all inclinations computed by 
RWIWD06. Fig. lc shows the histogram of the IRAC spectral indicies a for all 
observed point sources. The histograms of the IRAC a for YSOs in the Orion 
cloud and stars+YSOs in the IC348 region are shown for comparision. The 
colour and magnitude analysis shows that several point sources observed close 
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Figure 1. Histogram of the observed spectral indices (a) for sources associ- 
ated with candidate massive protostellar candidates. The dotted curve shows 
a similar histogram for the IRAC YSOs in Orion A cloud and the dashed 
curve shows the distribution of sources in IC348. 



to the HMPO targets display infrared excess emission in the IRAC bands that 
imitate the colours of YSO's in the Class I, II, III phases, and a good fraction 
of these are also luminous. Many point sources that are saturated in the IRAC 
bands may actually qualify as better IRCs. 37 saturated sources have been found 
to coincide well with the HMPO target positions in our sample. 

The colours and magnitudes analysis presented in Fig. 1 shows that: a) the 
point sources in the HMPO fields have unusually high a values (a >5), even 
higher than many of the deeply embedded Class I sources in Orion, showing 
that IRCs are indeed deeply embedded sources inside very dense cores and, b) a 
significant fraction of these IRCs are also luminous and represent 8-20 M Class 
I type objects when compared with accretion models. 

A total of 79 sources with AM>6 were classified as luminous IRCs to HMPO 
targets with photometry available in all four IRAC bands. Of these, 11, 27, 23 
and 25 are from the samples of Mol96, Sri02, Fon02, and Faun04 respectively and 
does not include the 37 saturated sources. These luminous IRCs are centred on 
the IRAS/mm peaks in the target fields with a nominal spread consistent with 
beam sizes and positional uncertainities. Some of these luminous IRCs can be 
reddened photospheres (particularly those with [5.8-8.0] <0. 4 colours) and some 
are probably evolved protostars; therefore, a complete SED analysis described 
in Sec. 4 was made. 



3. Infrared Nebulae 



The IRAC images of the HMPO targets reveal compact nebulae (10"-60" angular 
sizes) around several sources. The nebulae are found to be brightest in the 8 fim 
band and mostly invisible in the 3.6 (im IRAC band. We have used the 4.5, 
5.8 and 8.0 /im band images coded as blue, green, and red, respectively to 
generate composite false colour images for the observed nebulae. The analysis 
of the images are described by iKumar & Gravel (12007T ) and are also available 
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onlin^]]. The infrared nebulae repeatedly display cometary (e.g.: 118437-0216), 
core-halo (e.g: 118337-0743, 119403+2258), shell-like (e.g.: 119198+1423) and 
bipolar (e.g.:I18530+0215, 1192 13+1723) morphol ogies which are similar to the 
morphologies of UCHII regions (jKurtz et al.lll994l ) and bipolar outflows . Some 
images display only a single well-defined nebula, whereas others show a group 
of compact nebulae. A histogram of the projected sizes of these nebulae dis- 
plays a distribution in the range 0.1-1.0 pc with a mean value of 0.5 pc. These 
dimensions are similar to or smaller than the size of dense cores traced by the 
1.3 mm continuum maps dBeuther et alJl2002h . The 8.0 jim. band is dominated 



by PAH emission, which is known to be a tracer of radiation temperature. Ion- 
ising radiation from young massive stars, that may not yet be strong enough to 
produce a significant ionised region may therefore be traced by these infrared 
nebulae. Indeed a recent investigation has shown that the underlying structure 
of the ISM in such nebulae can possibly be inferred usin g the morphology o f 
the nebulae at various density regimes and ionising fluxes (jHeitsch et alJl2007h . 



Therefore, the morphology of the nebulae found here may well indicate the mor- 
phology in which ionising radiation is escapi ng from the unde r lying set up of 
physical structures close to the star. Recently IChurchwell et a l. (2006) used the 



GLIMPSE images to identify bubbles around OB stars in the Galaxy and argue 
that the smaller bubbles around several B stars are those produced by relatively 
soft radiation that fails to produce significant HII regions. The infrared neb- 
ulae around the HMPOs may well represent such bubbles or could be simple 
reflection nebulae due to an evolved generation of B stars. 



4. Radiative Transfer Modelling of the SEDs 

As described in Sec. 2, the high AM product sources are our best candidates 
of luminous IRCs to the HMPO targets. A bonafide list of 50 IRCs was made 
from the 79 high AM IRCs by using the constraint that the high AM product 
IRC was found within 2"of the observed mm/submm peak. We then extended 
the SEDs of this bonafide sample to the near-infrared and millimetre wavelength 
range by using the 2 MASS point source catalog, millimeter ob servations from 
Beuther et al. ( 20021 ) and t he submillimeter da ta from SCUBA ( Williams et al 



2004! ) and SIMBA cameras (jBeltran et al.l l2006). Data with the highest available 



angular resolution was used for the SEDs. In two cases, interferometric observa- 
tions in the mm and submm range were available that was utilised. The IRAC 
and 2MASS fluxes are restricted to their PSF sizes of less than 2"and typical 
mm/submm beam widths are ~6". The resulting SEDs were fi tted with an on- 
line S EP fitter tool developed and tested on low mass YSOs bv lRobitaille et al.l 
(120071 1. 



4.1. The SED fitting tool 

The SED fitting tool described by iRobitaille etaD (120071 ) is based on a grid of 



radiative transfer models described by RWIWD06. The models assume axisym- 
metric young stellar objects with a photosphere, disk, and infalling envelope. 



lr The colour images can be accessed at http://www.astro.up.pt/~nanda/hmpo/ 
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Figure 2. The best fit SED model for the source 1RAS18460-0307 and an 
8 /xm grey scale image of the same. The solid line shows the best fit model for 
an aperture of 15000AU, the dashed line represents disk emission, dot-dashed 
line the envelope component and the multiple-dot-dashed line the scattered 
component of the emission. Filled circles are data points and diamond symbols 
are upper limits. 



The upper limit of the mass range is 5OM0and the physics t hat governs the 
higher mass regime is based on accretion models described by IWhitney et al 



(|2004l 1. A 3re-calculated grid of model SEDs, computed using 14 variables are 



compared with the observed data. The models for which the difference between 
the x 2 value per data point of the fit and the best \ 2 value per data point is 
smaller than three are returned by the tool as the best fit models. Multiple 
models are returned by the tool for a given set of data points depending on 
the data quality, such as number of data points and their errors that define the 
SED. This results in a degeneracy which is typically 200 for the input SEDs 
we have used. In the presence of such high degeneracy, we adopt a method of 
constructing histograms for each physical parameter and finding the peak or by 
plotting a two parameter plot for all models and choosing the median of the 
distribution. For example to choose the best value of mass and time, we would 
plot a two parameter plot of mass vs time given by all the models, and if it 
results in a distribution centered on a particular mass and time value, those 
values are taken as the median representative va lues for mass and time of the 
source. These methods are explained in detail by I Grave & Kumar! (|2008ft . 



4.2. SED fitting results 

Figure. 2 shows the best fit model for the source IRAS 18460-0307 along with an 
8 /im image of the source. The solid line shows the best fit model for an aperture 
of 15000AU, the dashed line represents disk emission, dot-dashed line the enve- 
lope component and the multiple-dot-dashed line the scattered component of the 
emission. The median values of the physical parameters from model fits, such 
as mass, disk accretion rate, envelope accretion rate, stellar radii and age were 
tabulated for all the 50 sources. In those cases where higher spatial resolution 
data from the interferometers was available for the millimeter or submm bands, 
the degeneracy of the model fits was very low, demonstrating the importance of 
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such observations at longer wavelengths. Infact, a single model fit was obtained 
for the source IRAS18089-1732. 

The estimated physical parameters from the model fits span a range con- 
sistent with the known physics of massive stars. These parameters have a range 
of M = 6-45 M©, R = 10-100 R , t = 10 3 -10 5 yr, M disk = 0-0.9 M , and 
M env ~ 10 MqVt — lO~ 2 M0yr . These values are consistent with the input 
physics assumed in the SED model fitter. While the observed data is well mod- 
elled by assumed accretion scenario model, it is necessary to examine which of 
the estimated parameters are independent of the limitations of the model. In 
Fig. 3 we plot the accretion rate versus age of the source, representing the mass 
of the sources using symbols of different sizes. Circles represents disk accretion 
rate whereas triangles and squares represent envelope accretion rate. Squares 
are sources without disks. It can be seen from the figure that the most mas- 
sive stars (biggest symbols, Mgi3OM ) are all squares, meaning that they do 
not pocess any disks and are also among the youngest sources in the plot. The 
remaining symbols are more or less evenly spread around in time but not in 
M. Although the model grid consists of many models for massive stars with 
disks and envelopes together, the observed data best matches with massive star 
models without disks. Therefore, the result that "envelope accretion is the dom- 
inant factor in building massive stars" and the value M env ~ 10 -4 -10~ 2 Moyr -1 
is independent of any limitations that could be inherent in the model grid. 



5. Implications 

The SEDs of the IR counterparts to HMPO candidates can be well fitted by ac- 
cretion models involving an infalling envelope, disk, photosphere, bipolar cavity 
and a radiative equilibrium solution. In this scenario, the HMPO IRCs imitate 
the classical YSO phases known from low mass star formation. The envelope, 
rather than the disk, appears to be the major reservoir feeding the central engine 
in the most massive young stars. Therefore probing the inner few thousand AU 
structure of such source s are indispensable in understanding the details of mas- 
sive star formation (e.g. Grave Kumar 20071 ). 27 sources with IR counterparts 



also have some level of unresolved centimeter continuum emission which may rep- 
resent the ionised emission arising close to the central engines either in the form 
of expanding or infalling ionised regions and/or from winds. The presence of 
ionised gas together with an accretion set-up involving star/disk/envelope, with 
high accretion rates, suggest the possibility of the ionised gas in inward motions 
or forming a stable zone surrounding the central engine. Therefore it ap pears 
that t he th eoretical scenarios involving both ionised and molecular inflows (jKetol 
20021 . 120031 ) may well be the dominant mechanism producing the most massive 



stars. 



Reflection As evidenced by this meeting, there is an immediate need to better 
define the terminology such as cores, protostars and evolutionary stages in the 
massive star formation process. In the light of the results presented here, it 
appears that the envelopes around massive young stars are much larger than 
evidenced in the low mass star formation and plays an important role as mate- 
rial resorvoir to feed the central engine. The even larger flattenned structures/ 
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Figure 3. Accretion rate versus age plot from the best fit models of 50 
sources. Circles represents disk accretion rate whereas triangles and squares 
represent envelope accretion rate. Sizes of the symbols are proportional to 
the mass. 
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toroids observed in the millimeter regime, with sizes of 0.05 to O.lpc may repre- 
sent the basic units which collapse to form the massive stars. Thus, the definition 
of a Core may be better suited for such flattenned structures/toroids. We pro- 
pose an evolutionary sequence definition as follows: Stage I objects: millimeter 
bright without an IR point source (e.g: IRDC's and mm only cores); Stage II: 
IR bright (mm, IR, maybe weak cm continuum as well) with dominant envelope 
accretion, central engine is not directly contributing to the observed SED, high 
accretion rates and almost absent disks (sources from this study); Stage III: IR 
sources with lower accretion rates, near-infrared em ission lines , cent ral engine 
shows up in the near-infrared (e.g: sources from the Bik et ail (|2005l ) sample) 



Finally, we note that high angular observations of massive outflows needs to 
focus on investigating the ionised component and launching zones because the 
outflows from massive young stars may arise from the larger envelopes rather 
than disks according to the results of this study. 
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